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Colonization of a territory by a stochastic population under a strong Allee effect and 

a low immigration pressure 

Shay Be’er, Michael Assaf, and Baruch Meerson 
Racah Institute of Physics, Hebrew University of Jerusalem, Jerusalem 91904, Israel 

We study the dynamics of colonization of a territory by a stochastic population at low immigration 
pressure. We assume a sufficiently strong Allee effect that introduces, in deterministic theory, a large 
critical population size for colonization. At low immigration rates, the average pre-colonization 
population size is small thus invalidating the WKB approximation to the master equation. We 
circumvent this difficulty by deriving an exact zero-flux solution of the master equation and matching 
it with an approximate non-zero-flux solution of the pertinent Fokker-Planck equation in a small 
region around the critical population size. This procedure provides an accurate evaluation of the 
quasi-stationary probability distribution of population sizes in the pre-colonization state, and of the 
mean time to colonization, for a wide range of immigration rates. At sufficiently high immigration 
rates our results agree with WKB results obtained previously. At low immigration rates the results 
can be very different. 
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I. INTRODUCTION 

Any isolated population, which regulates itself via ran¬ 
dom births and deaths, is doomed to extinction [1-3]. 
Large and therefore long-lived stochastic populations ul¬ 
timately go extinct via a rare sequence of events when 
random population losses dominate over gains. This 
basic extinction scenario, unaccounted for by determin¬ 
istic theory, is at work in many situations in physics, 
chemistry, biology and other fields. One example from 
epidemiology is extinction of an endemic disease from 
a population when no new infectives arrive [1]. Prior 
to extinction, a large population resides in a long-lived 
quasi-stationary state, with a lifetime (the mean time to 
extinction) which is exponentially large in the average 
population size [3]. 

It has been long recognized that extinction is prevented 
by immigration: via either colonization of empty regions, 
or the “rescue effect” [4, 5]. Similarly, arrival of new in¬ 
fected individuals can restart the epidemics in a popula¬ 
tion which has recovered from an infection. Mathemati¬ 
cally, by introducing a constant immigration flux into a 
stochastic population model, one eliminates the absorb¬ 
ing state at zero population size and therefore prevents 
extinction. An important additional effect that many 
populations exhibit is the Allee effect, by which popula¬ 
tion biologists mean a group of effects causing a reduction 
in the per-capita growth rate at small population sizes [6] . 
In the language of deterministic theory, a strong Allee 
effect introduces a non-zero critical population size for 
establishment. If there is no immigration, and the initial 
population size is smaller than the critical size, the pop¬ 
ulation goes extinct quickly. If the initial population size 
is greater than the critical size, a long-lived state with 
a large population size appears. We will call this state 
the colonization state. In the presence of low immigra¬ 
tion pressure (by which ecologists mean small immigra¬ 
tion rate) the absorbing zero-population state gives way 


to the pre-colonization state: a state with a small pop¬ 
ulation size. As a result, the population can be either 
in the pre-colonization state, or the colonization state. 
The demographic noise (which, for large populations, is 
weak) causes rare switches between the two states. We 
will assume that the Allee effect is sufficiently strong so 
that this critical population size is large, see Sec. II for 
details. Here we evaluate the mean time to colonization 
(MTC), which we define as the mean switching time be¬ 
tween the pre-colonization and colonization states. 

Noise-induced switching between long-lived states is 
a classical paradigm of statistical physics going back to 
Kramers [7]. In the context of population dynamics, de- 
scribable as a continuous-time Markov process with a 
discrete space of states, accurate and useful general ex¬ 
pressions for the mean switching time have only become 
available recently. For single-population systems with 
single-step processes (that is, when there are only tran¬ 
sitions between a state with n individuals and a state 
with nil individuals), an exact analytical expression 
for the (properly defined) mean switching time can be 
obtained [8] by solving a recursive equation for the mean 
first passage time, see Sec. IV. This expression, however, 
is extremely cumbersome and not very informative. Fur¬ 
thermore, for multiple-step processes no exact solutions 
are available. These difficulties may explain the common 
practice, especially in the population biology literature 
[9], of using the so called “diffusion approximation”. In 
this approximation the mean switching time is evaluated 
from a Fokker-Planck equation that is derived from the 
original master equation via a truncated system-size ex¬ 
pansion [10]. Unfortunately, the diffusion approximation 
breaks down in the tails of the quasi-stationary distribu¬ 
tion. As shown in many studies [3, 11-15], this leads to 
errors in the mean switching time that are exponentially 
large in the population size, thus invalidating the whole 
calculation. 

A robust and efficient way of evaluating the mean 
switching time in large populations is provided by a dissi- 
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pative variant of WKB approximation [16] that employs 
the average population size as a large parameter and is 
applied directly to the original master equation [17]. In 
this way Dykman et al [18] calculated the effective en- 
tropic barrier that determines the mean switching time 
up to a pre-exponential factor. More recently, Meerson 
and Sasorov [19] used the WKB formalism to calculate, 
for a specific model, the mean switching time with ac¬ 
count of the pre-exponential factor. This calculation re¬ 
quires going to the sub-leading order of the WKB the¬ 
ory and also dealing with a vicinity of the unstable fixed 
point where WKB theory breaks down. The approach of 
Ref. [19] was extended to a general set of reactions by 
Escudero and Kamenev [20], see also Ref. [21]. 

The WKB formalism, however, assumes that the av¬ 
erage population size in each of the two bistable states 
is large. In the colonization problem this assumption 
breaks down when the immigration pressure is so small 
that the average population size in the pre-colonization 
state is 0(1) or less. Indeed, in some ecological systems 
the population dynamics represents a series of recur¬ 
rent extinctions and colonizations (the “rescue effect”) 
[4, 5, 22, 23]. In this case, and in other, milder, cases 
[9, 24, 25], the immigration rate is very small, and the 
situation calls for approximations that would respect 
the non-WKB character of the pre-colonization state. 
Here we develop an approximate method that yields 
the MTC and the long-lived quasi-stationary distribution 
(QSD) of population sizes in the situation when the pre¬ 
colonization state is of a non-WKB nature. The method 
holds for a broad range of immigration rates. For single- 
step processes it does not employ the WKB approxima¬ 
tion altogether, whereas for multiple-step processes the 
WKB approximation is only used for sufficiently large n 
where it is justified. 

An important element of our method is the zero-flux 
solution of the quasi-stationary master equation which 
can be found by recursion. For single-step processes, the 
resulting recursion solution is well known [8] , and it gives 
a very good approximation of the quasi-stationary distri¬ 
bution of the pre-colonization state for the population 
sizes n from n = 0 to a close vicinity of the Alice thresh¬ 
old (that is, of the unstable fixed point of deterministic 
theory). For multiple-step processes a recursion solution 
can often be obtained for sufficiently small n [21], and 
matched with a zero-flux WKB solution that remains 
valid until close to the Allee threshold. 

In the vicinity of the unstable fixed point, and at larger 
n , there is a finite probability flux toward larger n [19]. A 
proper non-zero-flux solution can be found by performing 
a boundary-layer analysis of the Fokker-Planck equation 
which can be derived from the original master equation 
and is valid in the vicinity of the Allee threshold [19]. By 
matching the zero-flux solution with the boundary-layer 
solution in their joint region of validity, we determine 
the probability flux and evaluate the MTC. The result¬ 
ing MTC exhibits an entropic barrier, so that the MTC is 
exponentially long. We show that the WKB approxima¬ 


tion remains remarkably accurate at quite low immigra¬ 
tion pressures, well beyond conservative estimates. As 
the immigration pressure decreases, the WKB prediction 
starts to fail. At not too small immigration rates, the 
correct entropic barrier still coincides with that obtained 
from WKB approximation. There is an important pre¬ 
factor, however, that strongly depends on the immigra¬ 
tion rate and, for very low immigration pressure, is very 
different from that predicted by the WKB approxima¬ 
tion. At still lower immigration rates even the entropic 
barrier is different from the WKB prediction. Our result 
for the MTC in a broad range of immigration rates is the 
central result of this work. 

For simplicity, we will present our method for a con¬ 
crete stochastic population model. Here is a plan of 
the remainder of the paper. The model is introduced 
in Sec. II. In Sec. Ill we present a derivation of the 
QSD of the pre-colonization state and obtain the MTC. 
In Sec. IV we compare our result for the MTC with the 
(very cumbersome) exact expression and with a WKB 
formula. The main results are summarized and discussed 
in Sec. V, while the Appendix contains a derivation of 
the Fokker-Planck equation and its approximate solution 
in the vicinity of the Allee threshold. 


II. MODEL 


We consider a stochastic population describable by 
a continuous-time and discrete-state Markov process. 
When only single-step processes are present, the master 
equation reads 


dP n 

^ — An —l-^n —1 A nPn T ^n-\-lPn-\-l l^nPn-i (1) 


where P„(f) is the probability of observing the popula¬ 
tion size n (n = 0,1,...) at time t, while A„ and fi n are 
the effective birth and death rates, respectively. The de¬ 
terministic rate equation, corresponding to the master 
equation (1), is 


,. — A n fi n . (2) 

at 

The specific model [26] we will be dealing with incorpo¬ 
rates an Allee effect, as modeled by Dennis [27], and a 
steady immigration flux, in a variant of the stochastic 
Verhulst model [21, 28]. In this model 


Bn 2 Bn 2 

A„ = r-|-—— and fj, n = n + ——, 

n + A K 


( 3 ) 


where time is rescaled so that the linear term in the 
death rate is equal to 1. The effective birth rate A„ ac¬ 
counts for immigration with n-independent rate r. The 
n-dependent part of \ n is proportional to n at large n 
and to n 2 at small n. Together with the linear in n part 
of the death rate this feature accounts for an Allee 
effect. The coefficient B > 1 is the reproduction rate. 
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The large parameters N and K control the Allee thresh¬ 
old and the carrying capacity of the colonization state, 
respectively. 

We will work in the parameter regime where Eq. (2) 
with rates (3) has three positive fixed points m < 
77-2 < 773 . The fixed points n\ and 773 are attracting. 
They correspond, in the deterministic theory, to the pre¬ 
colonization and colonization states, respectively. The 
fixed point 712 is repelling; it determines the Allee thresh¬ 
old. We will assume throughout the paper that the im¬ 
migration is weak, r -C N. In this regime the fixed point 
n\ can be obtained by neglecting the nonlinear terms in 
Eqs. (3) , while the other two fixed points can be obtained 
by neglecting the immigration: 


ni ~ r, 



1 

B 
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K 
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The deterministic equation (2) ignores demographic 
noise. The latter causes the population to switch ran¬ 
domly between the pre-colonization and colonization 
states. Figure (1) shows a typical realization of the 
stochastic dynamics of the system obtained in a Monte 
Carlo simulation employing the Gillespie algorithm [29] 
with the rates given in Eq. (3). The initial condition 
is such that the population finds itself, with probabil¬ 
ity close to 1, in the pre-colonization state. One can 
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FIG. 1. (Color online) Solid line: a typical Monte-Carlo real¬ 
ization of the random birth-death process n(t) with the birth 
and death rates given by Eq. (3). The parameters are B = 2, 
r = 2, N = 30 and K = 5000. The initial condition is n = 2. 
The dashed line shows the fixed point n = from Eq. (4). 
Inset: a blowup of n(t) in the pre-colonization state. The 
solid line depicts the simulation, the dashed line shows the 
fixed point n = ni from Eq. (4). 


see that the population dwells over a long time in the 
pre-colonization state. However, when a rare large fluc¬ 
tuation brings the population over the Allee threshold 
ri 2 , the population size flows almost deterministically to¬ 
wards the colonization state at n = 773 . Our task is to 
determine the QSD of the pre-colonization state and the 
MTC. 

For simplicity, we will assume that K is so large that 
the nonlinear term in the death rate /i n is negligible [30] . 
With the new rates, 

Bn? 

An = r H- —- and /x n = 77 , (5) 

77 + N 


the fixed point 77 = 773 moves to infinity, and the switch¬ 
ing problem is replaced by an effective problem of noise- 
driven population explosion. In this problem the popula¬ 
tion size, once it overcomes the Allee threshold, blows up 
in a finite time. This time scale is of deterministic nature 
and therefore relatively short [19]. The fixed points n\ 
and 772 become 


771 

N 

772 

N 


r 
N L 
1 


1 + 


B — 1 



and 


1 - B v+° 



( 6 ) 


respectively. 


III. SOLUTION 


A. Recursive solution 


When starting from a sub-threshold initial condition, 
the stochastic population relaxes, with a high probabil¬ 
ity, to the pre-colonization state around the stable fixed 
point ?7 1 of the deterministic theory. Although long-lived, 
this state is metastable, as there is a nonzero probability 
flux through the unstable fixed point 772 towards large 77 . 
At times much longer than the deterministic relaxation 
time (let us call it t r ), the pre-colonization probability 
distribution is described by the eigenvector of the master 
equation with the smallest positive eigenvalue 1/r, where 
r is the MTC [14, 19-21]: 

P n (t) ~ n n e~ t/T . (7) 

Here TT n is the QSD of the pre-colonization state. Plug¬ 
ging Eq. (7) into Eq. (1) and neglecting the exponentially 
small term — 7 r n /r on the left hand side, one arrives at a 
stationary difference equation for the QSD [14, 19-21]: 

An— l^n— 1 A n7T n T /in+l^n+l l^n^n = 0- (3) 

The disregard of the term —7 t u /t can only be justified 
if the immigration rate r is much larger than 1/t: a 
criterion that can be checked a posteriori. 
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Equation (8) is exactly soluble via recursion [8], and 
the zero-flux solution has the form 


'Kn 


-oil 


k -0 


Afc 

Mfc+1 


( 9 ) 


where ttq is determined from normalizing the total prob¬ 
ability to one [31]. For Afc and /j,k from Eq. (5) we obtain, 
with a help of “Mathematica”, 



. (10) 

where T(...) is the gamma function. This solution is 
valid for all n from zero to a close vicinity of the unsta¬ 
ble fixed point n 2 , see Refs. [19-21] and the next sec¬ 
tion. Employing the weak immigration assumption (that 
we have assumed in Sec. II), that is, r -C N, the pre¬ 
colonization QSD is sharply peaked around n\ — r. As a 
result, we can neglect, for the purpose of normalization, 
the second term, proportional to B, in the birth rate (5). 
The remaining simple immigration-death process is de¬ 
scribed by the Poisson distribution with mean r, and we 
find the normalization constant 7r 0 — e~ r . 

The distribution tail is described by the full expres¬ 
sion (10) that we will now simplify using the strong in¬ 
equalities n ~ N 1. Let us introduce the rescaled 
population size q = n/N that can be treated as a con¬ 
tinuous variable. We use the Stirling formula k! — 
y/2n k k+1 / 2 e~ k for the factorials of Eq. (10). For the 
squared absolute value of the gamma function in the nu¬ 
merator of Eq. (10) we can write |T(a + ib) | 2 = T(a + 
ib)T(a — ib). Using the Stirling formula in each of the 
multipliers, we obtain 


|T(o + ib) I 2 ~ 27 re _2a (o 2 + 7 , 2 )a-l/ 2 e - 2 &arcta„(b/a) _ 
After some algebra, all this yields 


r ( ( ?) — 


rA^/2ir(q + 1) 


BVW 

x e w [9 ln (^:)- ln (9+l)]-i[s+f-ln(JVg)] ; 


where 


A = 


T [ 1+ 2B +l 


Nr r 2 
~B ~ 4R2 


( 12 ) 


We will also need a more specialized asymptotic of 
Eq. (11) in a close vicinity of the unstable fixed point: 
]\r -1 / 2 -C (72 — q 1, where <72 = n 2 /N. Here it suffices 
to expand the logarithm of 7 r(q) in powers of q — q 2 up 
to second order: In [ 77 ( 9 )] = ao + ai(q 2 — q) + 02(92 — <z) 2 - 
As the coefficients oi and 02 are multiplied by a factor 


proportional to q — (72 "C 1, it suffices to calculate them 
only in the leading order in N ys> 1. The coefficient a 0 
demands a higher accuracy, and contributions of the or¬ 
der of O(N), 0(N 1 / 2 ) and 0(1) need to be kept. By 
doing so, and keeping terms up to 0(1) in the exponent, 
we can approximate Eq. (11) as 

\[2/ktA(B — 1 ) 

7T (q) ~ -;=- 

y/BN 

-Jv| ln ( B S_)-^.[ 2 B- 1 +ln(^ r i)] + i2_Jl_( 9 2-q) 2 | 
xe t > . (13) 


B. Boundary-layer solution 

The zero-flux approximation for the QSD, found in the 
previous section, is invalid close to the unstable fixed 
point 712 and at larger n, where the proper solution has 
a non-zero flux [19]. A major simplifying factor here 
is the validity of the Fokker-Planck approximation in 
a narrow boundary layer around ni'- \n — 712 ] "C 712 
or \q — < 72 1 "C 1 (to remind the reader, we assume 
B > 1). The Fokker-Planck approximation is valid 
here because the QSD varies sufficiently slowly with n: 
|(7r n +i — 7r n )/7r n | -C 1 [19-21]. The derivation of the 
Fokker-Planck equation and its solution in the boundary 
layer was presented elsewhere [19-21], For the reader’s 
convenience, we briefly reproduce these calculations in 
the Appendix. The boundary layer solution reads 


7T 


( BL )( q ) = (Q~<h 


■ e ^ erfc 


where 

2 _ 1 X(q 2 ) + n(q 2 ) 


Jc = 


Jc 


(14) 


(15) 


N (q 2 ) - /(q 2 y C (q 2 ) - n'(q 2 )' 

Here primes denote the derivatives with respect to the 
argument, A (q) = A n /N and n(q) = [i n /N. J c is the 
a priori unknown constant probability flux through the 
unstable fixed point n = n 2 . We will now determine it by 
matching the boundary layer solution (14) with the “bulk 
solution”, that is the zero-flux recursive solution (10), in 
their joint region of validity A 7-1 / 2 -C q 2 — q -C 1. The 
bulk solution is described in this region by the asymp¬ 
totic (13). Now we approximate the boundary layer so¬ 
lution (14) in this region. As l ~ A 7-1 / 2 , we can use the 
asymptotic erfc (—z) ~ 2 at z 1. Then, using Eqs. (5) 
and (6), we obtain after some algebra 

7 y BL \q) — V^ttBN J c e" 2® (<? 2 -«) _ (10) 


C. Quasistationary distribution and mean time to 
colonization 


Demanding that the expressions (13) and (16) coincide, 
we determine the probability flux J c : 

Jc ~ A &- i y e -N{H^)-M 2 *-^H^)]} (17) 

BN K ’ 
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FIG. 2. (Color online) The QSD 7 r„ as a function of n for 
B = 2, r = 2 and IV = 50 in a log scale. The QSD includes 
two overlapping asymptotics: the bulk asymptotic, given by 
Eq. (10) with 7ro = e~ r (solid line), and the boundary layer 
asymptotic, given by Eq. (14) with the flux J c from Eq. (17) 
(dashed line). 


where j(q) is the probability flux in the boundary layer. 
Taking the integral in Eq. (20), we obtain 


-N\j( q2 )-j( q .)}e-^. (21) 

We choose to satisfy 52 —<7- -A -1 / 2 , so that j(q~) — 

0. In its turn, j(q 2 ) ~ J c . Putting it all together, we 
obtain the MTC: 


B 


NJ C A(B - 1 )r 


0 iv{ln( -g^y ) + -g W [2B- l+ln( )]} 


( 22 ) 

This expression, valid in a broad range of immigration 
rates r, is the central result of this paper. It is instructive 
to consider different limits when this expression can be 
simplified. They are determined by the parameter 7 = 
sjNr/B that enters Eq. (12) for A. We obtain 


( 1 , 7 < 1 , 

\ |T(l + i 7 )| 2 , 7 = 0 ( 1 ), 

l 27776 ^ ln7_71 ’ 7 , 7>1. 


(23a) 

(23b) 

(23c) 


The QSD is now fully determined. The bulk of the QSD, 
for 0 < n < 7i2, is given by the zero-flux recursive solution 
Eq. (10) with 7 r 0 = e~ r . In the boundary layer |7z — 7t, 2 | -C 
n 2 the QSD is given by Eq. (14) with J c from Eq. (17) 
and A from Eq. (12). Figure (2) shows a plot of the QSD 
for a specific choice of parameters. 

We are now in a position to determine the MTC. Let 
us return to Eq. (1) for the time-dependent probability 
P n {t ) and sum it over n from n = 0 to n = 712 ■ As P„(i), 
at times t t r , is described by Eq. (7), the left hand 
side becomes 


The corresponding asymptotics 


T ~ 


B * n H- B^t) 


r{B - 1) 

|T(1 + i^)\ 2 B jvi n (- 
r(B - 1) 

27r 7 B j VIn f^B_ r ' ) _ 7rT 

r(B - 1) 


of the MTC are 
r 1/N, 

, r = 0(1/N), 


r > 1/N. 


(24a) 

(24b) 


(24c) 


E 


n —0 


dP n 

dt 


1 

T 



n —0 



(18) 


IV. COMPARISON WITH THE EXACT AND 
WKB RESULTS 


where we have used the fact that 7r n is normalized to 1, 
and the normalization is mostly contributed to by rela¬ 
tively small n’s. The summation over the right hand side 
of Eq. (1) can be split into two parts: 

ri 2 n— n 2 

£<-> = B->+£(•■•>. < 19 > 

71 = 0 71=0 71=71_ 


where n_ < n. 2 . Let us choose n_ so that it satisfies 
the double strong inequality V 1 / 2 -C n 2 — n_ N, or 
N -1 / 2 -C c /2 — q~ 1, where = ri-/N. Because of the 
inequality V 1 / 2 n 2 —the probability flux at any n < 

n_ is approximately zero, see the previous subsection. As 
a result, the first sum on the right hand side of Eq. (19) is 
zero. In its turn, the inequality 712 — n— N guarantees 
the applicability of the the Fokker-Planck approximation 
on the interval n_ < n < ri 2 - Therefore, the second term 
on the right in Eq. (19) can be approximated as 



dj(q) 

dq 


dq , 


( 20 ) 


A. Exact solution 

As we already mentioned, for single step processes the 
(properly defined) MTC can be found exactly from the 
backward master equation [8]. Although quite cumber¬ 
some, the exact solution is useful for our purposes, as it 
enables us to test the accuracy of our approximate result, 
Eq. (22). 

The exact derivation supposes that a single-step 
stochastic process with the birth rate A ra and death rate 
fi n is confined to the interval a < n < 6, where a and b 
are the reflecting and absorbing boundaries, respectively. 
The exact solution depends on the initial value n of the 
stochastic process. The mean time r(n) for the process 
to be absorbed at n = b obeys the exact equation [8] 

A n [r(n + 1) - r(n)] + fj, n [r(n - 1) - r(n)] = —1, (25) 
which should be solved with the boundary conditions 
r(a — 1) = r(a) and r(b) = 0. 


Ne~ t/T 
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In our case, A„ and p, n are taken from Eq. (3), while 
a = 0. A reasonable choice of b is the closest integer to 
the colonization fixed point n .3 = K(1 — 1 /B), whereas 
the initial number of individuals i is set to be the closest 
integer to r. The solution is [ 8 ] 

n 3 3 

r“ act (0 = E^)El/[ A ^)]> ( 26 ) 

j=i fc=0 


where 


m = ng- 


1 =1 


Although the products and sums in Eq. (26) can be 
brought to hypergeometric functions, it is more practical 
to evaluate Eq. (26) numerically. Figures 3 and 4 show 
comparisons of the exact result with that predicted by 
Eq. (22) at different r and A. As one can see, excellent 
agreement is observed for all relevant values of parame¬ 
ters. We also checked that the agreement is insensitive 
to the choice of the initial value n in the exact solution, 
as long as n < tt -2 and sufficiently far from n- 2 - 
We also compared Eq. (22) with results of extensive 
Monte Carlo simulations (not shown) and found excellent 
agreement. 


B. WKB approximation 

It is assumed in the existing formulations of the WKB 
theory that all relevant fixed points scale with the pop¬ 
ulation size [17-21]. In practice, the WKB approxima¬ 
tion is expected to hold in our colonization problem as 



r 


FIG. 3. (Color online) The ratio of the MTCs as a function 
of r for B = 2 and N = 50 in a log-log scale. Dashed line: 
T /r exact , where r is given by Eq. (22) and r exact j s given by 
Eq. (26). Solid line: t/t wkb , where r WKB i s given by the 
asymptotic limit r < Af of Eq. (23) of Ref. [20]. Dashed- 
dotted line: our theoretical prediction of t/t wkb = l /( 27 r 7 ) 
in the limit of 7 ^ 1, where t| 7<k i is given by Eq. (24a) and 
r WKB \ y< ^i is the 7 4C 1 asymptotic of Eq. (23) of Ref. [20]. 



FIG. 4. (Color online) The MTC as a function of A in a log 
scale for B = 5 and r = 1, 0.1, 0.01, and 0.001 in (a), (b), 
(c), and (d), respectively. Solid line: r from Eq. (22). Dashed 
line: r exact from Eq. (26). 


long as the pre-colonization fixed point corresponds to a 
sufficiently large population, even if r <C A. It is inter¬ 
esting to find out how large the pre-colonization popula¬ 
tion should be for the WKB theory to be accurate. We 
achieved this goal by comparing our approximate result 
for the MTC with that obtained via WKB approxima¬ 
tion. 

To calculate the MTC in the WKB approximation, 
t wkb , in the limit of r <C JV, we used Eq. (23) of Ref. 
[20] with the rates given by Eq. (5). We took the rates in 
the leading WKB order. As r dominates the birth rate 
in the vicinity of the pre-colonization fixed point, it has 
to be included in the leading WKB order. We simpli¬ 
fied the result by employing the smallness of the param¬ 
eter r/N. These calculations show that r n KB coincides 
with Eq. (24c) in the limit of r > 1/A. This inequal¬ 
ity is much weaker than the naively expected condition 
r 3 > 1 that would guarantee that the pre-colonization 
fixed point n\ ~ r is describable by a deterministic the¬ 
ory. It is surprising that the WKB theory remains accu¬ 
rate at much lower immigration pressures than one could 
have expected [32]. 

For r < 1/A, the WKB approximation breaks down. 
This is clearly seen in Fig. 3, where we compare our re¬ 
sult (22) for the MTC with the exact result and with 
the WKB result. In the limit of very low immigration 
pressure, r <C 1/A, there is a large factor missed by the 
WKB theory. This factor can be written as 


r WKB 


r<<1 ^ < 27 > 

7 or r that is exponentially small in the parameter A, the 
actor (27) describes an effective increase of the entropic 
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barrier to colonization, thus invalidating the WKB ap¬ 
proximation in its entirety [33]. 


V. SUMMARY AND DISCUSSION 
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We have investigated the dynamics of colonization of a 
territory by a stochastic population at low immigration 
pressure. Against all odds, and regardless of how small 
the immigration rate is, demographic noise eventually 
drives the population to the colonization state via a rare 
fluctuation that allows the population to overcome the 
Allee threshold. 

The specific model [26] we have dealt with incorporates 
an Allee effect, as modeled by Dennis [27], and a steady 
immigration flux, in a variant of the stochastic Verhulst 
model [21, 28]. We have determined the quasi-stationary 
distribution (QSD) of the population sizes and the mean 
time to colonization (MTC) in a broad range of immigra¬ 
tion pressures. In all parameter regions, our result for the 
MTC is in excellent agreement with the exact result and 
with Monte Carlo simulations. At moderate and high 
immigration rates our results agree with the previously 
found WKB results. At low immigration rates we obtain 
a large preexponential correction to the WKB result, due 
to the breakdown of the latter at low immigration rates. 
The correction factor becomes huge, and invalidates the 
WKB result completely, for very low immigration rates. 

The calculation method that we have presented here 
is free of uncontrolled assumptions and can be used for 
a broad class of stochastic population models that ex¬ 
hibit an Allee effect and low colonization pressure. It can 
also be extended to multi-step processes that, in general, 
do not admit exact solutions. In those cases the master 
equation can be linearized at small population sizes in 
order to determine the recursive solution there. The re¬ 
cursive solution can be matched, at n > 1 , with a bulk 
solution obtained by using the (leading and subleading 
order of the) WKB approximation [19-21]. The WKB 
approximation breaks down in the vicinity of the Allee 
threshold. There one has to match the WKB solution 
with the universal boundary-layer solution, obtained by 
solving the pertinent Fokker-Planck equation in a close 
vicinity of the Allee threshold. This double matching 
procedure is applicable for a broad class of multi-step 
processes involving arbitrary immigration rates, and it 
yields the QSD in the entire region of interest, and the 
MTC. 

In this work we have focused on colonization due to 
demographic stochasticity. It would be interesting to in¬ 
vestigate the interplay of demographic and environmental 
stochasticity [34-37] in the colonization under an Allee 
effect. It would be also interesting to study coloniza¬ 
tion under an Allee effect when new immigrants arrive in 
groups rather than separately [26, 38]. 


APPENDIX 


We briefly present here a derivation and solution of the 
quasi-stationary Fokker-Planck equation in the vicinity 
of the deterministic Allee threshold n = n 2 , see Sec. II. 
Let us define A( 9 ) = A n /N and fi(q) = /x„/A, where A„ 
and /j, n are the effective birth and death rates, respec¬ 
tively, given by Eq. (5). In terms of the rescaled variable 
q = n/N the quasi-stationary master equation ( 8 ) be¬ 
comes 


X(q - 1/N)n(q - 1/N) - X(q)n(q) 

+ m(<? + 1/N)n(q + 1/N) - = 0. (Al) 

Let us denote /+(<?) = X(q)Tr(q) and f-(q) = n(q)Tr(q). 
Taylor-expanding f± (9 =F 1 / N ) around q we find 

f±(q T 1/N) ~ f±(q) T f±(q)/N + /£(«)/(27V 2 ). (A2) 

Plugging Eq. (A2) into (Al) we arrive at the quasista¬ 
tionary Fokker-Planck equation 

lf + (q) - f-(q)} - ^[f+(q) + /-(<?)] = 0 , (A3) 

that can be written as d q J(q) = 0, where J is the proba¬ 
bility flux. Thus, the solution to this differential equation 
is a constant-flux solution. In order to proceed we notice 
that 7 r ( q) ~ Nn(q). Therefore, f + (q) is governed in the 
leading order by X(q)Tr (q ), while f_(q) is governed by 
n(q)TT ( 9 ). Integration over q yields, in the leading order 
of 1/N < 1, 

[•%) - Kq)Mq) - 7 ^( 9 ) + M<?)K (?) = Jo, (A4) 

where J c is the constant probability flux through the 
unstable fixed point. We can simplify Eq. (A4) in the 
boundary layer \q 2 — q\ <C 1 (see Sec. Ill B for the defini¬ 
tion of the boundary layer) by expanding the drift term 
up to the first order in q — q 2 and putting 9 = 92 in 
the diffusion term. The resulting equation accepts the 
universal form [19-21] 

{q ~ 92 ) 77 ( 9 ) - y 7 r '(9) = Jc, (A5) 


where i and J c are defined in Eq. (15) of the main text. 
Solving the linear first-order Eq. (A5), we obtain 


r = 


(12 gr 


( 9 ) = Ce & + 


Vnj c ( 12 -I) 2 

—e erf 


92 


(A 6 ) 






where C and J c are constants yet to be determined. C 
can be found by demanding that the boundary layer so¬ 
lution behaves properly at q > (72 [19]. Indeed, by con¬ 
sidering the asymptotic of the solution at q — 92 l we 


find that C = y/nJc/t in order to eliminate a rapid expo¬ 
nential growth. As a result, the boundary layer solution 
takes the form of Eq. (14) of the main text. 


[1] M. S. Bartlett, Stochastic Population Models in Ecology 
and Epidemiology (Wiley, New York, 1961). 

[2] R. M. Nisbet and W. S. C. Gurney, Modelling Fluctuating 
Populations (Wiley, New York, 1982). 

[3] O. Ovaskainen and B. Meerson, Trends in Ecology and 
Evolution 25 , 643 (2010). 

[4] J.H. Brown and A. Kodric-Brown, Ecology 58 , 445 
(1977). 

[5] I. Hanski, Metapopulation Ecology (Oxford University 
Press, Oxford, 1999). 

[ 6 ] P. A. Stephens, W. J. Sutherland, and R. P. Freckleton, 
Oikos 87 , 185 (1999); B. Dennis, ibid 96 , 3 (2002); F. 
Courchamp, J. Berec, and J. Gascoigne, Allee Effects in 
Ecology and Conservation (Oxford University Press, New 
York, 2008). 

[7] P. Hanggi, P. Talkner, and M. Borkovec, Rev. Mod. Phys. 
62 , 251 (1990). 

[ 8 ] C.W. Gardiner, Handbook of Stochastic Methods (Berlin, 
Springer, 2004). 

[9] A. Potapov and H. Rajakaruna, J. Theor. Biol. 337, 1 
(2013). 

[10] N. G. van Kampen, Stochastic Processes in Physics and 
Chemistry (North-Holland, Amsterdam, 2001). 

[11] B. Gaveau, M. Moreau, and J. Toth, Lett. Math. Phys. 
37, 285 (1996). 

[12] C. R. Doering, K. V. Sargsyan, and L. M. Sander, Mul¬ 
tiscale Model. Simul. 3, 283 (2005). 

[13] M. Assaf and B. Meerson, Phys. Rev. Lett. 97 , 200602 
(2006). 

[14] D. A. Kessler and N. M. Shnerb, J. Stat. Phys. 127, 861 
(2007). 

[15] M. Assaf and B. Meerson, Phys. Rev. E 75, 031122 
(2007). 

[16] C. M. Bender and S. A. Orszag, Advanced Mathemati¬ 
cal Methods for Scientists and Engineers (Springer, New 
York, 1999). 

[17] R. Kubo, K. Matsuo, and K. Kitahara, J. Stat. Phys. 9 , 
51 (1973). 

[18] M. I. Dykman, E. Mori, J. Ross, and P. M. Hunt, J. 
Chem. Phys. 100. 5735 (1994). 

[19] B. Meerson and P. V. Sasorov, Phys. Rev. E 78, 060103 
(2008). 


[20] C. Escudero and A. Kamenev, Phys. Rev. E 79 , 041149 
(2009). 

[21] M. Assaf and B. Meerson, Phys. Rev. E 81, 021116 

( 2010 ). 

[22] R.H. MacArthur and E. O. Wilson, Evolution 17 , 373 
(1963). 

[23] R.H. MacArthur and E. O. Wilson, The Theory of Is¬ 
land Biogeography (Princeton University Press, Prince¬ 
ton, 1967). 

[24] M. Vanhellemont, K. Verheyen, L. DeKeersmaeker, K. 
Vandekerkhove, and M. Hermy, Biol. Invas. 11, 1451 
(2009). 

[25] S. Dey and A. Joshi, Sci. Rep. 3, 1405 (2013). 

[26] B. Meerson and O. Ovaskainen, Phys. Rev. E 88 , 012124 
(2013). 

[27] B. Dennis, Natural Resource Modeling 3, 481 (1989). 

[28] I. Nasell, J. Theor. Biol. 211 , 11 (2001). 

[29] D.T. Gillespie, J. Phys. Chem. 81, 2340 (1977). 

[30] As we checked, the criterion for neglecting the nonlinear 
term in the death rate /x n is \[K N. 

[31] This normalization involves the summation of Eq. (9) 
over n in the region of 0 < n < ri 2 - 

[32] The technical reason for the good accuracy of the WKB 

asymptotic at r N -1 is the rapid convergence of the 

factor A from Eq. (12) to its leading-order Stirling ap¬ 
proximation at r A -1 . 

[33] To remind the reader, our formalism, based on the cal¬ 
culation of the QSD, demands that r be still much larger 
than the 1 /r. 

[34] A. Kamenev, B. Meerson, and B. Shklovskii, Phys. Rev. 
Lett. 101 , 268103 (2008). 

[35] E. Y. Levine and B. Meerson, Phys. Rev. E 87, 032127 
(2013). 

[36] M. Assaf, E. Roberts, Z. Luthey-Schulten, and N. Gold- 
enfeld, Phys. Rev. Lett. Ill, 058102 (2013). 

[37] M. Assaf, M. Mobilia and E. Roberts, Phys. Rev. Lett. 
Ill, 238101 (2013). 

[38] S. Be’er, M. Heller-Algazi and M. Assaf, Phys. Rev. E 
93, 052117 (2016). 



